import cv2
import numpy as np
import matplotlib.pyplot as plt
import math
import matplotlib
import os
def hilditch(img):
    # get shape
    H, W, C = img.shape

    # prepare out image
    out = np.zeros((H, W), dtype=np.int32)
    out[img[..., 0] > 0] = 1

    # inverse
    out = 1 - out

    while True:
        s1 = []
        s2 = []

        # step 1 ( rasta scan )
        for y in range(1, H - 1):
            for x in range(1, W - 1):

                # condition 1
                if out[y, x] > 0:
                    continue

                # condition 2
                f1 = 0
                if (out[y - 1, x + 1] - out[y - 1, x]) == 1:
                    f1 += 1
                if (out[y, x + 1] - out[y - 1, x + 1]) == 1:
                    f1 += 1
                if (out[y + 1, x + 1] - out[y, x + 1]) == 1:
                    f1 += 1
                if (out[y + 1, x] - out[y + 1, x + 1]) == 1:
                    f1 += 1
                if (out[y + 1, x - 1] - out[y + 1, x]) == 1:
                    f1 += 1
                if (out[y, x - 1] - out[y + 1, x - 1]) == 1:
                    f1 += 1
                if (out[y - 1, x - 1] - out[y, x - 1]) == 1:
                    f1 += 1
                if (out[y - 1, x] - out[y - 1, x - 1]) == 1:
                    f1 += 1

                if f1 != 1:
                    continue

                # condition 3
                f2 = np.sum(out[y - 1:y + 2, x - 1:x + 2])
                if f2 < 2 or f2 > 6:
                    continue

                # condition 4
                # x2 x4 x6
                if (out[y - 1, x] + out[y, x + 1] + out[y + 1, x]) < 1:
                    continue

                # condition 5
                # x4 x6 x8
                if (out[y, x + 1] + out[y + 1, x] + out[y, x - 1]) < 1:
                    continue

                s1.append([y, x])

        for v in s1:
            out[v[0], v[1]] = 1

        # step 2 ( rasta scan )
        for y in range(1, H - 1):
            for x in range(1, W - 1):

                # condition 1
                if out[y, x] > 0:
                    continue

                # condition 2
                f1 = 0
                if (out[y - 1, x + 1] - out[y - 1, x]) == 1:
                    f1 += 1
                if (out[y, x + 1] - out[y - 1, x + 1]) == 1:
                    f1 += 1
                if (out[y + 1, x + 1] - out[y, x + 1]) == 1:
                    f1 += 1
                if (out[y + 1, x] - out[y + 1, x + 1]) == 1:
                    f1 += 1
                if (out[y + 1, x - 1] - out[y + 1, x]) == 1:
                    f1 += 1
                if (out[y, x - 1] - out[y + 1, x - 1]) == 1:
                    f1 += 1
                if (out[y - 1, x - 1] - out[y, x - 1]) == 1:
                    f1 += 1
                if (out[y - 1, x] - out[y - 1, x - 1]) == 1:
                    f1 += 1

                if f1 != 1:
                    continue

                # condition 3
                f2 = np.sum(out[y - 1:y + 2, x - 1:x + 2])
                if f2 < 2 or f2 > 6:
                    continue

                # condition 4
                # x2 x4 x8
                if (out[y - 1, x] + out[y, x + 1] + out[y, x - 1]) < 1:
                    continue

                # condition 5
                # x2 x6 x8
                if (out[y - 1, x] + out[y + 1, x] + out[y, x - 1]) < 1:
                    continue

                s2.append([y, x])

        for v in s2:
            out[v[0], v[1]] = 1

        # if not any pixel is changed
        if len(s1) < 1 and len(s2) < 1:
            break

    out = 1 - out
    out = out.astype(np.uint8) * 255

    return out
def neighbours(x,y,image):
    """Return 8-neighbours of image point P1(x,y), in a clockwise order"""
    img = image
    x_1, y_1, x1, y1 = x-1, y-1, x+1, y+1
    return [ img[x_1][y], img[x_1][y1], img[x][y1], img[x1][y1], img[x1][y], img[x1][y_1], img[x][y_1], img[x_1][y_1] ]
def getSkeletonIntersection(skeleton):
    """ Given a skeletonised image, it will give the coordinates of the intersections of the skeleton.

    Keyword arguments:
    skeleton -- the skeletonised image to detect the intersections of

    Returns:
    List of 2-tuples (x,y) containing the intersection coordinates
    """
    # A biiiiiig list of valid intersections             2 3 4
    # These are in the format shown to the right         1 C 5
    #                                                    8 7 6
    validIntersection = [[0,1,0,1,0,0,1,0],[0,0,1,0,1,0,0,1],[1,0,0,1,0,1,0,0],
                         [0,1,0,0,1,0,1,0],[0,0,1,0,0,1,0,1],[1,0,0,1,0,0,1,0],
                         [0,1,0,0,1,0,0,1],[1,0,1,0,0,1,0,0],[0,1,0,0,0,1,0,1],
                         [0,1,0,1,0,0,0,1],[0,1,0,1,0,1,0,0],[0,0,0,1,0,1,0,1],
                         [1,0,1,0,0,0,1,0],[1,0,1,0,1,0,0,0],[0,0,1,0,1,0,1,0],
                         [1,0,0,0,1,0,1,0],[1,0,0,1,1,1,0,0],[0,0,1,0,0,1,1,1],
                         [1,1,0,0,1,0,0,1],[0,1,1,1,0,0,1,0],[1,0,1,1,0,0,1,0],
                         [1,0,1,0,0,1,1,0],[1,0,1,1,0,1,1,0],[0,1,1,0,1,0,1,1],
                         [1,1,0,1,1,0,1,0],[1,1,0,0,1,0,1,0],[0,1,1,0,1,0,1,0],
                         [0,0,1,0,1,0,1,1],[1,0,0,1,1,0,1,0],[1,0,1,0,1,1,0,1],
                         [1,0,1,0,1,1,0,0],[1,0,1,0,1,0,0,1],[0,1,0,0,1,0,1,1],
                         [0,1,1,0,1,0,0,1],[1,1,0,1,0,0,1,0],[0,1,0,1,1,0,1,0],
                         [0,0,1,0,1,1,0,1],[1,0,1,0,0,1,0,1],[1,0,0,1,0,1,1,0],
                         [1,0,1,1,0,1,0,0]];
    image = skeleton.copy()
    image = image/255
    intersections = list()
    for x in range(1,len(image)-1):
        for y in range(1,len(image[x])-1):
            # If we have a white pixel
            if image[x][y] == 1:
                neighbour = neighbours(x,y,image)
                valid = True;
                if neighbour in validIntersection:
                    intersections.append((y,x))
    # Filter intersections to make sure we don't count them twice or ones that are very close together
    for point1 in intersections:
        for point2 in intersections:
            if (((point1[0] - point2[0])**2 + (point1[1] - point2[1])**2) < 10**2) and (point1 != point2):
                intersections.remove(point2)
    # Remove duplicates
    intersections = list(set(intersections))
    return intersections
def degree(x1, y1, x2, y2):
    angle = 0.0
    dx = x2 - x1
    dy = y2 - y1
    if x2 == x1:
        angle = math.pi / 2.0
        if y2 == y1:
            angle = 0.0
        elif y2 < y1:
            angle = 3.0 * math.pi / 2.0
    elif x2 > x1 and y2 > y1:
        angle = math.atan(dx / dy)
    elif x2 > x1 and y2 < y1:
        angle = math.pi / 2 + math.atan(-dy / dx)
    elif x2 < x1 and y2 < y1:
        angle = math.pi + math.atan(dx / dy)
    elif x2 < x1 and y2 > y1:
        angle = 3.0 * math.pi / 2.0 + math.atan(dy / -dx)
    return angle * 180 / math.pi
def mkdir(path,name):
    # 去除首位空格
    path=path.strip()
    # 去除尾部 \ 符号
    path=path.rstrip("/")
    list2 = path.split('/')
    strr=""
    for i in range(len(list2)-1):
        if i==0:
            strr=list2[0]
        else:
            strr=strr+"/"+list2[i]
    s=strr+"/"+name
    #判断路径是否存在
    isExists=os.path.exists(s)
    # 判断结果
    if not isExists:
        # 如果不存在则创建目录,创建目录操作函数
        '''
        os.mkdir(path)与os.makedirs(path)的区别是,当父目录不存在的时候os.mkdir(path)不会创建，os.makedirs(path)则会创建父目录
        '''
        #此处路径最好使用utf-8解码，否则在磁盘中可能会出现乱码的情况
        os.makedirs(s.encode("utf-8").decode('utf-8'))
        print(s+' 创建成功')
        #return True
    else:
        # 如果目录存在则不创建，并提示目录已存在
        print(s+' 目录已存在')
        #return False
"""
CNFD函数接口：
调用方式：CNFD(file_pathname,filename)
输入：参数分别是file_pathname（文件路径名）以及filename（文件名）  eg:CNFD("D:/mywork/Arbitrary/gt","FJ103.png")
输出：返回结果是连通域的个数，同时会在文件路径名下生成一个Color文件夹，保存不同连通域的分类效果
"""
def CNFD(file_pathname,filename):
    img = cv2.imread(file_pathname + '/' + filename)
    img1 = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(img1, connectivity=8)
    # print(filename)
    # print('主连通域个数', num_labels - 1)
    CNFD= num_labels - 1
    #print('Nerve Fibre Density(CNFD)= ', CNFD)  # 400x400 µm2=0.16mm2
    output = np.zeros((img1.shape[0], img1.shape[1], 3), np.uint8)
    for i in range(1, num_labels):
        mask = labels == i
        output[:, :, 0][mask] = np.random.randint(0, 255)
        output[:, :, 1][mask] = np.random.randint(0, 255)
        output[:, :, 2][mask] = np.random.randint(0, 255)
    mkdir(file_pathname, "Color")
    cv2.imwrite("./Color/" + filename, output)
    return CNFD
"""
CNBD函数接口：
调用方式：CNBD(file_pathname,filename)
输入：参数分别是file_pathname（文件路径名）以及filename（文件名）  eg:CNFD("D:/mywork/Arbitrary/gt","FJ103.png")
输出：返回分支点的个数，同时会在文件路径名下生成一个Refinement文件夹保存血管细化的结果，一个Branch_point文件夹保存血管分支点标记的结果
"""
def CNBD(file_pathname,filename):
    refinement = cv2.imread(file_pathname + '/' + filename).astype(np.float32)
    XHout = hilditch(refinement)
    mkdir(file_pathname, "Refinement")
    cv2.imwrite("./Refinement/" + filename, XHout)
    img = cv2.imread(file_pathname+'/'+filename)
    img11 = cv2.imread("./Refinement/" + filename, 0)
    img111 = getSkeletonIntersection(img11)
    #print("所有的分支点数", len(img111))
    CNBD=len(img111)
    #print("Nerve Branch Density(CNBD)", len(img111) / 0.16)
    for item in img111:
        data = np.array(item)
        cv2.circle(img, (data[0], data[1]), 2, (0, 0, 255), 3)
    mkdir(file_pathname, "Branch_point")
    cv2.imwrite("./Branch_point/" + filename, img)
    fenzhi = cv2.imread("./gt/" + filename, 0)
    for item in img111:
        data = np.array(item)
        cv2.circle(fenzhi, (data[0], data[1]), 3, (0, 0, 255), 6)
    mkdir(file_pathname, "Segmentation")
    cv2.imwrite("./Segmentation/" + filename, fenzhi)
    return CNBD
"""
CTBD函数接口：
调用方式：CTBD(file_pathname,filename)
输入：参数分别是file_pathname（文件路径名）以及filename（文件名）  eg:CNFD("D:/mywork/Arbitrary/gt","FJ103.png")
输出：返回分支数量
"""
def CTBD(file_pathname,filename):
    img = cv2.imread(file_pathname+'/'+filename)
    img1 = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(img1, connectivity=8)
    img11 = cv2.imread("./Refinement/" + filename, 0)
    img111 = getSkeletonIntersection(img11)
    imgsu = cv2.imread("./Segmentation/" + filename)
    imgss = cv2.cvtColor(imgsu, cv2.COLOR_BGR2GRAY)
    num__labels, labels, stats, centroids = cv2.connectedComponentsWithStats(imgss, connectivity=8)
    NUMBER=num__labels - num_labels - len(img111)
    CTBD=(num__labels - num_labels - len(img111)) / 0.16
    #print("所有的分支数", num__labels - num_labels - len(img111))
    #print("Nerve Fibre Total Branch Density (CTBD)", (num__labels - num_labels - len(img111)) / 0.16)
    return NUMBER
"""
CNFL函数接口：
调用方式：CNFL(file_pathname,filename)
输入：参数分别是file_pathname（文件路径名）以及filename（文件名）  eg:CNFD("D:/mywork/Arbitrary/gt","FJ103.png")
输出：返回所有血管的长度，同时会在文件路径名下生成一个Degree文件夹，保存血管角度的直方图统计，直方图以20度为一个间隔
"""
def CNFL(file_pathname,filename):
    img = cv2.imread(file_pathname + '/' + filename)
    imgsu = cv2.imread("./Segmentation/" + filename)
    imgss = cv2.cvtColor(imgsu, cv2.COLOR_BGR2GRAY)
    ret, thresh = cv2.threshold(imgss, 230, 255, cv2.THRESH_BINARY_INV)
    contours, hierarchy = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    sum = 0
    matrix = [0] * 180
    for c in contours:
        x, y, w, h = cv2.boundingRect(c)  # 计算点集最外面的矩形边界
        cv2.rectangle(imgss, (x, y), (x + w, y + h), (0, 255, 0), 2)
        (x, y), radius = cv2.minEnclosingCircle(c)
        center = (int(x), int(y))
        radius = int(radius)
        sum = sum + radius * 2
        cv2.circle(imgss, center, radius, (0, 0, 255), 2)
        # 对每个轮廓点求最小外接矩形
        rect = cv2.minAreaRect(c)
        # cv2.boxPoints可以将轮廓点转换为四个角点坐标
        box = cv2.boxPoints(rect)
        # 这一步可以保证四个角点坐标为顺时针
        startidx = box.sum(axis=1).argmin()
        box = np.roll(box, 4 - startidx, 0)
        tangle = degree(box[0][0], box[0][1], box[1][0], box[1][1])
        # print(int(tangle))
        matrix[int(tangle)] = matrix[int(tangle)] + 1
        rect_points = np.int0(box)
        cv2.drawContours(img, [rect_points], 0, (255, 255, 255), 2)
    sum = sum * 1.0417 * 0.001
    #print("各段总长", sum)
    #print("Nerve Fibre Length(CNFL)", sum / 0.16)
    return sum
    count = 0
    i = 0
    listt = [0] * 180
    lii = [0] * 10
    while (count < 180):
        summ = matrix[count:count + 20]
        numm = summ[0] + summ[1] + summ[2] + summ[3] + summ[4] + summ[5] + summ[6] + summ[7] + summ[8] + summ[9] + summ[
            10] + summ[11] + summ[12] + summ[13] + summ[14] + summ[15] + summ[15] + summ[16] + summ[17] + summ[18] + \
               summ[19]
        listt[count] = numm
        count = count + 20
        lii[i] = numm
        i = i + 1
    print(' ')
    plt.title('图中血管的角度统计')  # 添加标题
    matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
    x = range(0, 10)
    plt.xticks(x, ("0", "20", "40", "60", "80", "100", "120", "140", "160", "180"), color='blue')
    plt.bar(x, height=lii, label='角度统计', width=0.5, facecolor='#9999ff', edgecolor='white')
    plt.xlabel('角度')
    plt.ylabel('数量')
    mkdir(file_pathname, "Degree")
    plt.savefig('./Degree/' + filename)
    plt.close()





